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Abstract. Wave packets provide a well established and versatile tool for studying 
time-dependent effects in molecular physics. Here, we demonstrate the application of 
wave packets to mesoscopic nanodevices at low temperatures. The electronic transport 
in the devices is expressed in terms of scattering and transmission coefficients, which are 
efficiently obtained by solving an initial value problem (I VP) using the time-dependent 
Schrodinger equation. The formulation as an IVP makes non-trivial device topologies 
accessible and by tuning the wave packet parameters one can extract the scattering 
properties for a large range of energies. 



1. Introduction. 

Low-dimensional mesoscopic devices are systems in which electronic movement is 
confined along a two-dimensional layer within the semiconductor pQ. Molecular-beam 
epitaxy allows one to produce extremely clean systems where the free path length of 
electrons at very low temperature reaches from the fim up to the mm scale before 
inelastic collisions occur and coherence is lost [2] . The large coherence length requires to 
treat the system on a quantum-mechanical level since interference and diffraction effects 
occur. The electron densities in the electron layer are on the order of 10 15 m -2 , with a 
positively charged background of approximately the same density located at a distance 
of 30 — 100 nm from the electronic layer. The electrons are occupying the available states 
up to the Fermi energy Ep, which can be adjusted by external gates. Experimentally, 
the properties of the system are probed by alloyed metallic contacts which reach down 
to the electronic layer. At low temperatures the current through the system is measured 
as a function of the applied voltages at the contacts. In principle, the electric circuit is 
closed through a battery, where charges are separated and after travelling through the 
wires get injected in the semiconductor [3]. A quantum-mechanical description of the 
complete circuit cannot be achieved. Instead one implements an open system, where the 
contacts inject the electrons coming from an external reservoir, which can be viewed as 
an environment. The determination of the correct boundary conditions at the contacts 
is a considerable challenge and presents a still unsolved problem for interacting particles. 
Theoretically, transport in mesoscopic systems at low temperatures and at low currents 
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is often described within non-interacting models. The justification for this drastic 
simplification lies in the idea that only electrons with energies close to the Fermi energy 
are participating in the net transport and that Pauli blocking suppresses scattering 
events between electrons below and above the Fermi energy completely. Considerable 
effort has to be made to correctly determine the effective potential landscape for the 
electrons at the Fermi energy, since the positive background charges are largely screened 
by the electrons within the device. The effect of screening has to be calculated for a 
specific geometry in the presence of metallic gates and randomly distributed donor 
charges [TJ H]. The resulting mean- field approach to mesoscopic transport has worked 
reasonably well for a wide variety of systems, but interactions effects are known to 
play an important role in small devices and in magnetic fields, where charging effects 
(Coulomb blockade) and correlations (fractional quantum Hall effect) occur. 

For predicting the current through the system in response to a small change of 
voltage at one of the contacts, quantitative calculations rely on knowing the transmission 
(or scattering) probability from the source contact to the drain contact. The irregular 
device geometry necessitates the use of numerical methods to evaluate the scattering 
cross-sections. Scattering processes can be described in a stationary (time-independent) 
way using eigenstates of the asymptotic part of the system or using a time-dependent 
approach based on the propagation of wave packets. Here, we will discuss the time- 
dependent wave packet method, which allows us in principle to handle both, non- 
interacting and interacting systems. The wave packet techniques used in atomic and 
molecular physics have to be changed in several ways in order to be applicable for 
nanodevices [5]. We will first give a detailed discussion of the non- interacting theory 
and in the last section comment on interaction effects. 

2. Device layout and the asymptotic region 

The proper definition of a scattering matrix requires to introduce asymptotic regions 
where the scattering potential is absent and exact eigenstates can be constructed. In 
the two-dimensional systems under consideration, the contacts inject electrons into 
waveguide like structures which guide the electrons to the scattering region. If we 
consider the contacts as emitting incoherently electrons, we obtain the emitted electron 
current by convolution of the product of the local density of states (LDOS) at the 
point of injection with the group velocity and the occupation probability given by the 
Fermi-Dirac distribution. 

2.1. Parabolic waveguide in a homogeneous magnetic field 

An instructive example consists of an electron released into a parabolic potential 
V(y) = \ muj yV 2 along the ^-direction. In addition a homogeneous magnetic field B is 
present along the z-axis together with a strongly confining potential along the z-directon, 
which quantizes the system along this axis and reduces the effective dimensionality of 
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the subsystem under consideration to two dimensions. The Hamiltonian for the two- 
dimensional motion reads in the Landau gauge A = (— By, 0, 0) [6]: 



H = 2^ ( ^ + qBy)2 ~ 2^* + 2 muj y y2 



(1) 



Here, m denotes an effective electronic mass (in AlGaAs/GaAs layers m = 
0.067 m e ), which results from taking into account the periodic crystal structure of 
the semiconductor within the k • p expansion p] . It is useful to introduce the cyclotron 
frequency ujq = ff and to express u y in terms of another frequency Q 2 = u 2 + uJq. The 
normalized eigenfunctions and eigenenergies of the Hamiltonian ([!]) are: 
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We express the LDOS in terms of a sum and an integral over the squared eigenfunctions 
weighted with the eigenenergies: 
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For each quantum number n, the energy integrated LDOS becomes 
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No spatial dependence of N n remains after the energy integration. The density of 
electrons per quantized level (the degeneracy) differs from the quantization in a purely 
magnetic field (or orthogonal electric and magnetic fields) by the factor Q 2 /oJq [7]. 
We proceed to calculate the quantum mechanical current emitted with fixed energy 
from the point r [8J. The drift current j(r;E) is given by the product of the group 
velocity d k E nk /h times the local density of states [9]. We obtain two counterpropagating 
components of the current: 
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The total current through a closed surface around r is given by the sum j(r;E) = 
\j Xj+ (r; E)\ + \j x _(r; E)\. At energies %Q,(n + 1/2) the product of the singular LDOS 
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Figure 1. Current flow in a wave guide coupled to two reservoirs with Fermi energies 
4 0) and4 L) . 



with the group velocity results in a finite current. The energy integral over the current 
weighted with the occupation probability gives 
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where f(E, E F , T) = l/(exp[(J5 — E F )/(ksT)} + 1) denotes the Fermi-Dirac distribution. 
We obtain the global current by an integration over the local current given in eq. (jSJ): 

pL poo 

J{E F ) = dx dy J(r; E F ) 



Jo 

2Le 00 



OO 
OO 



2nh 



/OO 
dE f(E, E F , T)Q(E - hn(n + 1/2)). 



(10) 



From eq. (JTUl) . we calculate the current through a parabolic waveguide which is coupled 
at x = to electrons in a reservoir with Fermi energy E F [0,y] 
x = L to a second reservoir with Fermi energy £7jr[W, y] = E{ 



J F 



Ep and at 



(see fig. [T]). For 

simplicity, we set the temperature to zero. The net current J net along x — ... L is 
the sum of the counterpropagating currents from both reservoirs (each reservoir has 
an extension of AL along x): J+{Ef) = J(eP)/2 and J-(E { F L) ) = J(E { F L) )/2. For 
m(M + 1/2) < E^ < EP < hVL(M + 3/2), M = 0, 1,2, . . . 

j m - J{Ef)) : J{Ep) = e(M : 1} k> - 4 L) ) . (id 

net 2AL 2nh V F ' v ; 

where M + 1 is the number of occupied modes below the Fermi energy. The net current 
depends only on the difference of the Fermi energies along the waveguide. Eq. (fTT!) 
can be obtained without explicitly referring to the LDOS and the local drift current 
(see P, eq. (12.9)), but the apparent detour using the LDOS makes the connection 
with the source approach to quantum transport [8] more transparent. Since the squared 
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are normalized to l/(27i~), the spatial integration just gives a 



factor AL/ (2tt), while the partial derivative of E n ^ x cancels the inverse partial derivative 
coming from the energy integration of the DOS: 
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Figure 2. Two-terminal resistivity l/cr ne t as function of the magnetic field. The left 



reservoir is at E^p 



16.5 meV, the right one at Ep = 17.5 meV 



Eq. ({TO]) is converted to a conductivity by dividing J net by the difference in Fermi 
energies: 

ff "^ e ira^ (M + 1) ^ (13) 

t r 

The resulting conductivity is quantized in units of e 2 /h, as long as the voltage difference 
between the two reservoirs does not exceed the mode separation HQ. The resulting 
resistivity p = 1 / cT not is shown in fig. [2] for a fixed Fermi energy as function of the 
magnetic field. The origin for the quantization in steps of e 2 /h is not the quantized 
number of available states per Landau level ((6]), but rather the effective reduction to a 
longitudinal one-dimensional channel. 



2.2. Scattering and multiple channels 

The conductance quantization in eq. ( IT3|) holds only in the absence of scattering within 
the waveguide. In the presence of scattering due to deformations and impurities in 
the waveguide, the conductivity will get reduced. Numerical methods are required 
to calculate the transmission and conductance of such a system. The attachment of 
additional asymptotic channels (terminals) to the waveguide shown in fig. |3] requires to 
extend the formalism to handle the multi-terminal case. This generalization was done 
by Biittiker [10J and leads to the following expression for the current from channel i 
e r°° 

I i = T dE E {UnanjiE^ifiE^^-fiE^^T)), (14) 

where ti ni j nj denotes the transmission amplitude for scattering from the transverse mode 
Uj in arm j into the mode rii in arm i. 



3. From wave packets to transmissions 

The Landauer-Buttiker formalism relates the elements of the scattering matrices with 
the currents and voltages in the system. The system consists of a central region 
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Figure 3. Geometry of a four terminal device with arms 1 and 4 and main channels 
2 and 3. The simulation area covers a width of 0.7 /im and a height of 0.9 /xm. (a) 
The sender wave packet sits in the lower arm and the receiver wave packet in the left 
arm in order to track the time-dependent correlation function, (b) Flux lines in the 
asymptotic regions are indicated, through which the time-dependent probability flux is 
recorded, (c) Complex absorbing potentials (CAPs) at the end of the channels mimic 
a semi-infinite extension and are used to record the flux leaving the simulation region. 

of irregular shape, which is connected to semi-infinite channels, which are assumed 
to be free of imperfections which cause scattering. The semi-infinite channels allow 
one to introduce a set of channel eigenstates and to study transmission from a well- 
defined asymptotic state. A similar concept is behind Wigner's and Eisenbud's R- 
matrix approach, where an artificial boundary is introduced at which the asymptotic 
regions and the interaction/scattering region get connected. Within the system the 
probability current is conserved and all equations can be derived by considering the 
probability flux through a closed surface. The probability flux describes stationary [TT] 
as well as time-dependent processes [12J. The time-dependent approach transforms 
the scattering problem into an initial value problem (IVP), which circumvents the 
explicit construction of eigenstates. The adventage of the IVP formulation become most 
apparent if transmission coefficients for a wide energy range are required (to study how 
the conductivity changes with Fermi energy, non-zero temperatures, or small voltages). 
In practice we construct a wave packet of finite extent by forming a superposition of 
plane waves along the waveguide with a specific transversal mode of the waveguide. 

4. Propagation methods 

An analytic time-evolution of the wave packet in the scattering region is only known for 
very few problems [13J. Thus in general numerical methods are required to solve the 
time-dependent Schrodinger equation for this part of the system. Roughly the methods 
can be divided in two classes, (i) methods based on the Trotter expansion of the time- 
evolution operator and (ii) polynomial expansions of the time-evolution operator. The 
time evolution typically depends on a kinetic part T = —% 2 V 2 / (2m) and a potential part 
V(r) which do not commute. In the Trotter approach [H] the completete time interval 
is broken into N small time-steps At and one expands the time-evolution operator into 
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products of exponential functions for the kinetic and the potential part 

U(t, t + NAt) = e-^+^'/ft « (e-iVAt/2H e -iTAt/n e -iVAt/2ny (15) 

The commutation error of potential and kinetic energy in this expansion is of the 
order At 3 and can be neglected for small enough chosen time-steps. Usually the initial 
wave function ip(r,t) is given in position space and consequently the potential operator 
exp(— iVAt/2h) is diagonal in this representation. In order to apply the kinetic part 
exp(— iTAt/h) a Fourier transform J 7 is used to express the wave function in momentum 
space. Thus the total time evolution of the wave packet is given by 

^(r,t + NAt)^e-' lVAt/2h 

This scheme is easily adapted to various systems and gives converged results provided 
that the time-step At is small enough [5]. Another possibility, especially suited for 
long time scales, is the direct expansion of the time-evolution operator in Chebyshev 
polynomials C n [15] 

M 

ij(r,t + NAt) = J2^-Sn, )Jn(AENAt/{2h))C n (-iH Iloim )'iP(r,t), (17) 

n=0 ' ' 

■l n 

where J n denotes the Bessel function of order n. The order of the required polynomials 
M is increased until Jm{AEN At/ (2Ti)) drops below the desired precision. Since the 
Chebyshev polynomials form a perfect expansion of the exponential function on the 
interval [—1,1], the propagation converges in the energy range AE if the Hamiltonian 
is rescaled as H norra = 2T-L/AE. The states P n are obtained by the recursion relations 
P = ip(r,t), Pi = H norm ip(r,t), and P n = —2iH norm P n _i + P n -2- The commutator error 
is absent in the polynomial expansion and long propagation times can be achieved. 

In an open system the scattering region is coupled to asymptotic channels. One way 
to take into account the semi-infinite extent of the channels is to introduce a complex 
absorbing potential (CAP) iU(r) which is added to the potential V(r). CAPs are widely 
used in quantum chemistry [16]. The complex- valued potential U(r) is chosen to be zero 
in the scattering region and gradually grows in the channels towards the border of the 
simulation region (see fig. E^c)) in order to avoid spurious reflections. As a result the 
norm of the wave function decreases during the propagation. The complex term iU(r) 
violates the self-adjointness of the Hamiltonian and leads to complex eigenvalues. Since 
the high-order Chebyshev expansion does not converge away from the real axis Faber or 
Newton polynomials have to be used for long-term propagations in open systems [153 • 

Tracking the probability flux 

The transmission probabilities through the mesoscopic device is extracted from the time 
evolution of wave packets by recording the probability flux in the asymptotic regions. 
There are several possible ways to obtain the flux [17], (i) cross- correlation functions, 
(ii) fluxlines, or (iii) using complex absorber. The flux is recorded either time-resolved 
or energy resolved. Each approach has specific (dis)adventages and all of them can be 
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used in conjunction. In the cross correlation approach we define a sender wave packet 
describing the initial situation. It is located in channel i populating transversal mode 
Hi and is composed of purely incoming momenta directed towards the scattering region. 
The destination, channel j with transversal mode rij, is represented by an outgoing 
receiver state. We record the overlap between forward propagated sender wave packet 
"0send( r , t) and stationary receiver wave packet ip rec (r) 

C(t) = (VwlVWd(i))- (18) 

The transmission amplitude is then given by the energy representation of the cross 
correlation function C(t) 

W*0 = ~. (F\ r ( 19 ) 

The factors r] send and fi rec are correction terms which account for the longitudinal shape of 
the wave packets and ensure that t ini j nj (E) does not depend on wave packet parameters. 
A second way consists in calculating the probability flux through a line in the asymptotic 
channels. The flux operator is given by 

F = —(p8(x-x Q ) + 5(x-x )p), (20) 
2 m 

where x denotes the longitudinal direction of the channel and xq is the postion of the 
fluxline. The expectation value (ip\F\ifj) = J j(x ,y)dy gives the integrated flux along 
the fluxline. The transmission propability is obtained by 

\tin ijnj (E)\ 2 = (^mhnfP^l^E)), (21) 

where 

Kn, = J W,nSE')}^- ni (E')\dE' (22) 

projects onto outgoing states in channel i with transversal mode n,. Since the fluxline 
is located in the asymptotic channel we can relate the scattering eigenstate \ip~ n .(E)) 
to channel eigenstates ip a ,n,± = Xa,n,±k x e ±lkxXa . The scattering eigenstate 

\^U(E)} = r ^ s Ut)e iEt/n dt (23) 

is extracted from the propagation history of the sender wave packet. The third method 
works by calculating the reduction of the norm due to presence of the CAP within each 
channel and to obtain the flux from the change of the norm between successively applied 
absorption events in each channel. 



5. Application example: Hall cross 

As an example of the formalism, we discuss the transmission through the device shown 
in fig. |3l which is a four terminal device and can be used to probe the Hall effect. For a 
numerical evaluation, several practical constraints have to be fulfilled: the wave packet 
must fit on the numerical grid in configuration and in momentum space. In order to 
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Figure 4. Transform of the time-dependent cross-correlation to the energy-resolved 
transmission for the Hall cross in fig. [3] at a magnetic field of B = 0.4 T. (a) Time- 
dependent cross-correlation function for the transmission of lead 2. mode into arm 
1, mode 1. (b) Transmission amplitude, obtained by the Fourier transform of (a), (c) 
Current / through the system as function of Fermi energy for a voltage of 0.1 mV and 
temperature 1 K. The dashed line corresponds to the step-wise changes of the current 
in the harmonic waveguide, eq. (|12l) . (d) Resistivity as function of the Fermi energy. 

use the fast Fourier transform (FFT) it is desirable to use the same number of grid 
points in both representations [H]. The discretization of the position and momentum 
space gives a lower limit for the spatial extension of the used wave packets, since one 
has to ensure that each packet is sampled sufficiently. And last, numerically we cannot 
propagate the wave packets infinitely long so that we have to stop the time evolution 
after a certain time. The propagation time is chosen long enough that most parts of 
the wave packet already traveled through the device and got absorbed in the CAPs at 
the end of the channels. In fig. H^a) we illustrate the detected overlap between receiver 
and sender wave packets. The sender was chosen to represent the transversal ground 
state in lead 2 and the receiver populates the first exited mode in arm 1. During 
time evolution the sender packet propagates through the Hall cross, scatters in different 
channels, and finally gets absorbed. Hence the cross correlation function decays with 
time. After 120 ps the absolute value of the cross correlation remains below 5 x 10 -5 
and we stop the propagation. From eq. ( fl9l) . the transmission propability is obtained 
by performing the Fourier transform of the cross correlation function which is done 
efficiently with a FFT. We get the transmissions for discrete energy values where the 
spacing is determined by the total propagation time. For the present simulation we got 
116 energy points per meV. Fig. H^b) shows the propability for inter mode scattering 
from lead 2 with mode to arm 1 in mode 1. 

In contrast to quantum chemistry where the wave packet initially represents an 
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eigenstate of a certain potential energy surface, wave packets are not associated with a 
direct physical meaning in their application for mesoscopic transport. Rather they are 
used as tools to compute transmission probabilities for a scattering potential and we can 
match the longitudinal shape to the needs at hand. It is convenient to use Gaussian wave 
packets whose parameters are adjusted to represent a certain energy range. A single 
wave packet run gives then the transmission amplitudes for an energy range which is 
sufficiently represented by the used wave packet. To increase the numerical efficiency we 
set the Gaussian width as small as possible to get an extended energy representation. 
A further enhancement of efficiency is achieved by using several receiver states within 
a single propagation run. Here, we put five receiver states in each lead populating the 
transversal modes up to n = 4. Thus we obtain about 40,000 transmission amplitudes 
corresponding to a typical computation time of 0.2 s per amplitude on a standard CPU. 
Current and voltages in the four terminal Hall cross are evaluated by inserting the 
computed transmission amplitudes in eq. (fll]) . We apply a bias voltage of 0.1 mV 
which drops symmetrically between contact 2 and contact 3 around the Fermi energy 
Ep. The temperature is set to 1 K corresponding to 4/cgT = 0.35 meV. Contacts 1 and 
4 act as perfect voltage probes forcing the currents l\ and I4 to vanish. This gives a 
nonlinear system of equations, whose solution determines the chemical potentials /ii and 
/x 4 and hence the voltage drop V# = (fJ-i—^/e- In fig.|U(c) we show the current between 
lead 2 and lead 3 in dependency on the Fermi energy. The conductance quantization of 
the harmonically confined channel gives rise to a step like behavior of the current. Note 
that the current remains below the conductance quantization of a harmonic waveguide 
(dashed line) since electron transmission in the side arms as well as reflections diminish 
the conductance through the Hall cross. Furthermore the non-zero temperature smears 
out structures on the energy scale of 4k B T and especially any abrupt jumps at mode 
openings get broadened. Fig. @Jd) shows the Hall resistance Rh = Vh/I% which goes 
down with increasing number of open modes. Since single mode effects superpose each 
other, Rh = Vn/h reflects an averaged behavior and gets more regular with increasing 
Fermi energies. 

6. Open questions in mesoscopic physics 

The time-dependent approach to mesoscopic physics can also be generalized to materials 
with more complicated band-structures. In general one has to propagate wave packets 
with additional spin or pseudo-spin degrees of freedom, like the two-dimensional time- 
dependent Dirac equation for graphene [IB] . 

The previous calculation assumed an effective non-interacting description of the 
electrons in the nanodevice. The inclusion of interactions beyond the mean-field level 
is a complicated task, but is required in order to explain many observed phenomena. 
As an example we mention the quantum Hall effect, where the existence of the Hall 
potential is a direct consequence of interaction effects [HI [20]. In principle the wave 
packet formalism can be extended to incorporate interactions between electrons, but this 
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generalization requires to propagate a properly (anti-) symmetrized product of single- 
particle orbitals. A more serious conceptional problem arises due to the injection of 
electrons into the device, which requires to change the device many-body electronic wave 
function. Most current approaches switch off the interactions away from the scattering 
region. However, this procedure may not be compatible with boundary conditions at 
the contacts. Self-consistent schemes are under development but require further work 
in order to explain the experimental observations [211 HI 120] . 
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